//!
//!       Staggered Mesh for u-vel and v-vel               Notation
//!      
//!       0   1       2       3       4   5             
//!       |   |       |       |       |   |            o Center of volumes    
//!
//!   5-  #---^--->---^--->---^--->---^---#  -4        > u-velocity 
//!       |       |       |       |       |   
//!   4-  >   o   >   o   >   o   >   o   >            ^ v-velocity 
//!       |       |       |       |       |            
//!       ^---^---+---^---+---^---+---^---^  -3        + nodes
//!       |       |       |       |       |         
//!   3-  >   o   >   o   >   o   >   o   >            # corners, here >, ^
//!       |       |       |       |       |              and o are defined
//!       ^---^---+---^---+---^---+---^---^  -2 
//!       |       |       |       |       |            
//!   2-  >   o   >   o   >   o   >   o   >            
//!       |       |       |       |       |            
//!       ^---^---+---^---+---^---+---^---^  -1 
//!       |       |       |       |       |                 
//!   1-  >   o   >   o   >   o   >   o   >                 
//!       |       |       |       |       |                 
//!   0-  #---^--->---^--->---^--->---^---#  -0 
//!
//!      |       |       |       |       |
//!      0       1       2       3       4                  
//!
//!
//!         Indexation of u (>), v (^) and other variables (o).
//!                                                  
//!                      v(i,j) 
//!                  |     n     |               
//!                --+-----^-----+--             
//!                  |           |               
//!                  |           |               
//!     u(i-1,j) = w >     o     > e = u(i,j)
//!                  |   (i,j)   |
//!                  |           |
//!                --+-----^-----+--
//!                  |     s     |
//!                     v(i,j-1)
//!

namespace Tuna {

template<class Tprec, int Dim>
inline bool CDS_Hay<Tprec, Dim>::calcCoefficients1D() 
{
    prec_t G_dx = Gamma / dx;
    prec_t dx_dt = dx / dt;
    prec_t ce, cep, cem, cw, cwp, cwm, CE, CW;
    aE = 0.0; aW = 0.0; aP = 0.0; sp = 0.0;

    for (int i =  bi; i <= ei; ++i) 
      {
	CE = ce = u(i ) ;
	CW = cw = u(i-1) ;	 
	cem = cep = 0;
	cwm = cwp = 0;

	if ( ce > 0 ){
	  CE = 0;
	  cep = ce * 0.5 * (-phi_0(i) + phi_0(i+1));
	} else {
	  cem = ce * 0.5 * (phi_0(i) - phi_0(i+1));
	} 

	if ( cw > 0 ){
	  cwp = cw * 0.5 * (-phi_0(i-1) + phi_0(i));
	} else {
	  CW = 0.0;
	  cwm = cw * 0.5 * (phi_0(i-1) - phi_0(i));
	} 

	aE (i) = G_dx - CE;
	aW (i) = G_dx + CW;
	aP (i) = aE (i) + aW (i) + dx_dt + (ce - cw);
	sp (i) = phi_0(i) * dx_dt - (cep + cem - cwp - cwm);
      }
    applyBoundaryConditions1D();
    return 0;
}

template<class Tprec, int Dim>
inline bool CDS_Hay<Tprec, Dim>::calcCoefficients2D() 
{
    prec_t Gdy_dx = Gamma * dy / dx;
    prec_t Gdx_dy = Gamma * dx / dy;
    prec_t dxy_dt = dx * dy / dt;
    prec_t ce, cep, cem, cw, cwp, cwm, CE, CW;
    prec_t cn, cnp, cnm, cs, csp, csm, CN, CS;
    aE = 0.0; aW = 0.0; aN = 0.0; aS = 0.0; aP = 0.0; 
    sp = 0.0;

    for (int i =  bi; i <= ei; ++i)
      for (int j = bj; j <= ej; ++j)
	{
	  CE = ce = u(i  , j) * dy;
	  CW = cw = u(i-1, j) * dy;
	  CN = cn = v(i, j  ) * dx;
	  CS = cs = v(i, j-1) * dx;
	  cem = cep = 0;
	  cwm = cwp = 0;
	  cnm = cnp = 0;
	  csm = csp = 0;
	  
	  if ( ce > 0 ){
	    CE = 0;
	    cep = ce * 0.5 * (-phi_0(i,j) + phi_0(i+1,j));
	  } else {
	    cem = ce * 0.5 * (phi_0(i,j) - phi_0(i+1,j));
	  } 
	  
	  if ( cw > 0 ){
	    cwp = cw * 0.5 * (-phi_0(i-1,j) + phi_0(i,j));
	  } else {
	    CW = 0.0;
	    cwm = cw * 0.5 * (phi_0(i-1,j) - phi_0(i,j));
	  } 

	  if ( cn > 0 ){
	    CN = 0;
	    cnp = cn * 0.5 * (-phi_0(i,j) + phi_0(i,j+1));
	  } else {
	    cnm = cn * 0.5 * (phi_0(i,j) - phi_0(i,j+1));
	  } 
	  
	  if ( cs > 0 ){
	    csp = cs * 0.5 * (-phi_0(i,j-1) + phi_0(i,j));
	  } else {
	    CS = 0.0;
	    csm = cs * 0.5 * (phi_0(i,j-1) - phi_0(i,j));
	  } 
	  	    
	  aE (i,j) = Gdy_dx - CE ;
	  aW (i,j) = Gdy_dx + CW ;
	  aN (i,j) = Gdx_dy - CN ;
	  aS (i,j) = Gdx_dy + CS ;
	  aP (i,j) = aE (i,j) + aW (i,j) + aN (i,j) + aS (i,j) + dxy_dt
	    + (ce - cw) + (cn - cs);
	  sp (i,j) = phi_0(i,j) * dxy_dt
	    - (cep + cem - cwp - cwm + cnp + cnm - csp - csm);	    
	}
    applyBoundaryConditions2D();
    return 0;
}


//!
//!---------------------  3D  ---------------------
//!
template<class Tprec, int Dim>
inline bool CDS_Hay<Tprec, Dim>::calcCoefficients3D() 
{
    prec_t dyz = dy * dz, dyz_dx = Gamma * dyz / dx;
    prec_t dxz = dx * dz, dxz_dy = Gamma * dxz / dy;
    prec_t dxy = dx * dy, dxy_dz = Gamma * dxy / dz;
    prec_t dxyz_dt = dx * dy * dz / dt;
    prec_t ce, cep, cem, cw, cwp, cwm, CE, CW;
    prec_t cn, cnp, cnm, cs, csp, csm, CN, CS;
    prec_t cf, cfp, cfm, cb, cbp, cbm, CF, CB;
    aE = 0.0; aW = 0.0; aN = 0.0; aS = 0.0; aF = 0.0; aB = 0.0; aP = 0.0; 
    sp = 0.0;

    for (int k = bk; k <= ek; ++k)
      for (int i =  bi; i <= ei; ++i)
	for (int j = bj; j <= ej; ++j)
	  {
	    CE = ce = u(i  , j, k) * dyz;
	    CW = cw = u(i-1, j, k) * dyz;
	    CN = cn = v(i, j  , k) * dxz;
	    CS = cs = v(i, j-1, k) * dxz;
	    CF = cf = w(i, j, k  ) * dxy;
	    CB = cb = w(i, j, k-1) * dxy;
	    cem = cep = 0;
	    cwm = cwp = 0;
	    cnm = cnp = 0;
	    csm = csp = 0;
	    cfm = cfp = 0;
	    cbm = cbp = 0;
	  
	    if ( ce > 0 ){
	      CE = 0;
	      cep = ce * 0.5 * (-phi_0(i,j,k) + phi_0(i+1,j,k));
	    } else {
	      cem = ce * 0.5 * (phi_0(i,j,k) - phi_0(i+1,j,k));
	    } 
	  
	    if ( cw > 0 ){
	      cwp = cw * 0.5 * (-phi_0(i-1,j,k) + phi_0(i,j,k));
	    } else {
	      CW = 0.0;
	      cwm = cw * 0.5 * (phi_0(i-1,j,k) - phi_0(i,j,k));
	    } 	    
	    
	    if ( cn > 0 ){
	      CN = 0;
	      cnp = cn * 0.5 * (-phi_0(i,j,k) + phi_0(i,j+1,k));
	    } else {
	      cnm = cn * 0.5 * (phi_0(i,j,k) - phi_0(i,j+1,k));
	    } 
	    
	    if ( cs > 0 ){
	      csp = cs * 0.5 * (-phi_0(i,j-1,k) + phi_0(i,j,k));
	    } else {
	      CS = 0.0;
	      csm = cs * 0.5 * (phi_0(i,j-1,k) - phi_0(i,j,k));
	    } 

	    if ( cf > 0 ){
	      CF = 0;
	      cfp = cf * 0.5 * (-phi_0(i,j,k) + phi_0(i,j,k+1));
	    } else {
	      cfm = cf * 0.5 * (phi_0(i,j,k) - phi_0(i,j,k+1));
	    } 
	    
	    if ( cb > 0 ){
	      cbp = cb * 0.5 * (-phi_0(i,j,k-1) + phi_0(i,j,k));
	    } else {
	      CB = 0.0;
	      cbm = cb * 0.5 * (phi_0(i,j,k-1) - phi_0(i,j,k));
	    } 
	
	    aE (i,j,k) = dyz_dx - CE;
	    aW (i,j,k) = dyz_dx + CW;
	    aN (i,j,k) = dxz_dy - CN;
	    aS (i,j,k) = dxz_dy + CS;
	    aF (i,j,k) = dxy_dz - CF;
	    aB (i,j,k) = dxy_dz + CB;
	    aP (i,j,k) = aE (i,j,k) + aW (i,j,k)  + aN (i,j,k) + aS (i,j,k)
	      + aF (i,j,k) + aB (i,j,k) + dxyz_dt
	      + (ce - cw) + (cn - cs) + (cf - cb);
	    sp (i,j,k) = phi_0(i,j,k) * dxyz_dt 
	      - (cep + cem - cwp - cwm + cnp + cnm - csp - csm + cfp + cfm - cbp - cbm); 
	  }
    applyBoundaryConditions3D();   
    return 0;
}
  
} // Tuna namespace














